What’s inside atl object

Now, let’s take a look the atl object in which we read the GTFA feed from MARTA. This object is a list. In it, names(atl) shows that there are 9 data.frames. Notice that there are many tables displayed in Table 1 above, but only some of them are in atl. Notice in Table 1 that only some are required; others are optional.

typeof(atl)

## [1] "list"

names(atl)

## [1] "agency"         "calendar"       "calendar_dates" "routes"        
## [5] "shapes"         "stop_times"     "stops"          "trips"         
## [9] "."

print(head(atl))

## $agency
## # A tibble: 1 × 7
##   agency_id agency_name                  agenc…¹ agenc…² agenc…³ agenc…⁴ agenc…⁵
##   <chr>     <chr>                        <chr>   <chr>   <chr>   <chr>   <chr>  
## 1 MARTA     Metropolitan Atlanta Rapid … http:/… Americ… en      (404)8… custse…
## # … with abbreviated variable names ¹​agency_url, ²​agency_timezone,
## #   ³​agency_lang, ⁴​agency_phone, ⁵​agency_email
## 
## $calendar
## # A tibble: 4 × 10
##   service_id monday tuesday wednesday thursday friday saturday sunday start_date
##   <chr>       <int>   <int>     <int>    <int>  <int>    <int>  <int> <date>    
## 1 2               0       0         0        0      0        0      0 2021-08-14
## 2 3               0       0         0        0      0        1      0 2021-08-14
## 3 4               0       0         0        0      0        0      1 2021-08-14
## 4 5               1       1         1        1      1        0      0 2021-08-14
## # … with 1 more variable: end_date <date>
## # ℹ Use `colnames()` to see all variable names
## 
## $calendar_dates
## # A tibble: 3 × 3
##   service_id date       exception_type
##   <chr>      <date>              <int>
## 1 5          2021-09-06              2
## 2 5          2021-11-25              2
## 3 5          2021-11-26              2
## 
## $routes
## # A tibble: 118 × 8
##    route_id route_short_name route_lon…¹ route…² route…³ route…⁴ route…⁵ route…⁶
##    <chr>    <chr>            <chr>       <chr>     <int> <chr>   <chr>   <chr>  
##  1 15664    1                Marietta B… ""            3 https:… FF00FF  ""     
##  2 15665    2                Ponce de L… ""            3 https:… 008000  ""     
##  3 15666    3                Martin Lut… ""            3 https:… FF8000  ""     
##  4 15667    4                Moreland A… ""            3 https:… FF00FF  ""     
##  5 15668    5                Piedmont R… ""            3 https:… 00FFFF  ""     
##  6 15669    6                Clifton Ro… ""            3 https:… 008080  ""     
##  7 15670    8                North Drui… ""            3 https:… 008000  ""     
##  8 15671    9                Boulevard … ""            3 https:… 00FFFF  ""     
##  9 15672    12               Howell Mil… ""            3 https:… FF00FF  ""     
## 10 15673    14               14th Stree… ""            3 https:… 00FF00  ""     
## # … with 108 more rows, and abbreviated variable names ¹​route_long_name,
## #   ²​route_desc, ³​route_type, ⁴​route_url, ⁵​route_color, ⁶​route_text_color
## # ℹ Use `print(n = ...)` to see more rows
## 
## $shapes
## # A tibble: 370,776 × 4
##    shape_id shape_pt_lat shape_pt_lon shape_pt_sequence
##    <chr>           <dbl>        <dbl>             <int>
##  1 119402           33.8        -84.6                 1
##  2 119402           33.8        -84.6                 2
##  3 119402           33.8        -84.6                 3
##  4 119402           33.8        -84.6                 4
##  5 119402           33.8        -84.6                 5
##  6 119402           33.8        -84.6                 6
##  7 119402           33.8        -84.6                 7
##  8 119402           33.8        -84.6                 8
##  9 119402           33.8        -84.6                 9
## 10 119402           33.8        -84.6                10
## # … with 370,766 more rows
## # ℹ Use `print(n = ...)` to see more rows
## 
## $stop_times
## # A tibble: 1,033,990 × 5
##    trip_id arrival_time departure_time stop_id stop_sequence
##    <chr>   <time>       <time>         <chr>           <int>
##  1 8044645 14:35:00     14:35:00       903320              1
##  2 8044645 14:36:02     14:36:02       903448              2
##  3 8044645 14:36:51     14:36:51       901144              3
##  4 8044645 14:37:43     14:37:43       904219              4
##  5 8044645 14:38:19     14:38:19       903850              5
##  6 8044645 14:38:39     14:38:39       903664              6
##  7 8044645 14:39:39     14:39:39       213268              7
##  8 8044645 14:40:12     14:40:12       213101              8
##  9 8044645 14:41:35     14:41:35       212926              9
## 10 8044645 14:41:53     14:41:53       211886             10
## # … with 1,033,980 more rows
## # ℹ Use `print(n = ...)` to see more rows

Calculate travel times from one station to all other stations

The tidytransit package offers a very convenient function that calculates the shortest travel times from a stop to all other stops. Try ?tidytransit::travel_times to see what this function does and what arguments it takes. Notice that only two out of nine arguments are required, and the rest are optional. The required arguments are ‘filtered_stop_times’ and ‘stop_name.’ If you do not provide the optional arguments, the function will use the pre-populated values for the calculation.

But if you run the function right away, there will be an error. The error is caused by the fact that there are stations that share names but are far from each other. There may be more than one stations in different parts of the city that share their names. If this was the case, the duplicate names would not be errors in the data. However, if in reality there are no such stations, the duplicate names would be actual flaws in your data.

The travel_times() function calculates travel time based on station names. This function is designed to throw an error if there are stops with the same name but are more than 300 meters apart. We need to eliminate duplicate names through the following steps:

1. For each unique station names in stops table, calculate how many duplicates there are and, if there are duplicates, calculate distances to the duplicate stops. stop_group_distances() function can do this.
2. Get names of the stops that have duplicates far from each other. If stops with duplicate names are close (e.g., within 200 meters), they may actually be same stations.
3. For all stops with duplicates > 200-meter distances, change their station names by appending (1), (2), … (n) where n is the number of duplicates. For example, if there are two “MIDTOWN STATION”s, they will each be “MIDTOWN STATION (1)” and “MIDTOWN STATION (2)”.
4. One of the requirements for travel_times() function is that atl object contains a transfer table, which is absent when freshly imported from the hard drive. This is NOT always the case – the transfer table is an optional table in GTFS definition. Some entities include transfer table by default, some don’t. MARTA is among those who don’t include transfer table. We can use gtfs_transfer_table() function from gtfsrouter package to create this table and include it to the atl object.
5. Finally, one of the two required arguments, ‘filtered_stop_times’, takes an output from a function called tidytransit::filter_stop_times(). This function takes two required arguments: (1) our atl object and (2) a date for which we’d like calculate travel time. You can also choose to specify (3) minimum departure time and (4) maximum arrival time.

Finally, you are ready to calculate travel times!

# Step 1:
stop_dist <- stop_group_distances(atlsf$stops, by='stop_name') %>% 

# Step 2:  
  filter(dist_max > 200)

# Step 3
atl$stops <- atl$stops %>% 
  group_by(stop_name) %>% 
  mutate(stop_name = case_when(stop_name %in% stop_dist$stop_name ~ paste0(stop_name, " (", seq(1,n()), ")"),
                               TRUE ~ stop_name)) 

# Step 4
atl$transfers <- gtfsrouter::gtfs_transfer_table(atl, 
                                     d_limit = 200, 
                                     min_transfer_time = 120)
# Step 5
am_stop_time <- filter_stop_times(gtfs_obj = atl, 
                                  extract_date = "2021-08-14",
                                  min_departure_time = 3600*7, # input unit is in second. So 3500*7 is 7AM
                                  max_arrival_time = 3600*10)  # similarly, this is 10AM 

# travel_times
trvt <- travel_times(filtered_stop_times = am_stop_time,
                     stop_name = "MIDTOWN STATION",
                     return_coords = TRUE)

trvt_pal <- colorQuantile(palette = "Reds", domain = trvt$travel_time)

leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addCircles(data = trvt %>% st_as_sf(coords = c("to_stop_lon", "to_stop_lat"), crs = 4326),
             fillColor = ~trvt_pal(travel_time),
             stroke = F,
             radius = 300,
             fillOpacity = 0.7,
             popup = paste0("Travel Time from <br> <strong> MIDTOWN: ", round(trvt$travel_time/60, 2), " minutes <strong>") %>% 
               lapply(htmltools::HTML))

On average, you can see that it takes 50.1 minutes from the Midtown Station to other stations. On the map, you can see that it takes longer time to get to stops that are further away from the Midtown Station. Let’s do this calculation for other stops in different parts of the city to see if there are a particular corners of the city that are burdened with longer travel time on average. To do this, we will do the following steps:

1. Randomly sample a few stations from each Census Tract.
2. TBD

# Average time to get to other stations
print(str_c("On average, it takes ", 
            round(mean(trvt$travel_time)/60,1), 
            " minutes to travel from the Midtown Station to all other stations."))

## [1] "On average, it takes 50.1 minutes to travel from the Midtown Station to all other stations."

Calculate area of service by SES

GTFS from MARTA contains information about both rails and buses. The ‘routes’ table in atl has a column named ‘route_type,’ which contains integer values denoting types of transit service. To see what each integer means, you will need to look at the link to the Reference provided above (or see below for a screen capture from the Reference).

Reference from Google

‘routes’ table contains route_types, which is needed for separating rail transit. ‘shapes’ table contains the shapefile needed for geo-operations such as buffer, intersections, etc. Because ‘routes’ and ‘shapes’ table do not contain a common key, they need to be joined through ‘trips’ as intermediate table.

# Filter out rail transit
route_trip <- atl$trips %>% dplyr::left_join(atl$routes, by = "route_id")

unique_shape <- route_trip %>% 
  dplyr::group_by(route_id) %>% 
  dplyr::slice(1) %>% 
  dplyr::ungroup()

route.shape <- atlsf$shape %>% inner_join(unique_shape, by = "shape_id")

# Route type is not really intuitive - let's fix that
route.shape <- route.shape %>% 
  dplyr::mutate(route_type = dplyr::case_when(
    route_type == "0" ~ 'Tram, Streetcar',
    route_type == "1" ~ 'Subway, Metro',
    route_type == "2" ~ 'Rail',
    route_type == "3" ~ 'Bus'
  ))

pal <- leaflet::colorFactor(c("red", "orange", "pink"), domain = route.shape$route_type)

route.shape %>% 
  leaflet::leaflet(data = .) %>% 
    leaflet::addProviderTiles(providers$CartoDB.DarkMatter) %>% 
    leaflet::addPolylines(color = ~pal(route_type), 
                 weight = 3,
                 opacity = 0.9,
                 popup = paste0("Route type: ", route.shape$route_type))

Until Version 1.0, the sf package used ‘equirectangular projection’ for many operations that involve, for example, intersects, intersection, union, nearest_point, st_join, etc. It did so even when the given data is in geographical coordinates.

From Version 1.0, sf package uses S2 spherical geometry for spatial operations, based on S2 geometry library written by Google. For more information on what all of these means, you can read this post.

# If we do not convert CRS to a projected one, the buffer may generate a ragged boundaries. 
# https://r-spatial.github.io/sf/articles/sf7.html#buffers-1
# https://r-spatial.org/r/2020/06/17/s2.html#sf-10-goodbye-flat-earth-welcome-s2-spherical-geometry
sf::sf_use_s2(TRUE)
MARTA_buffer <- route.shape %>% 
  sf::st_transform(crs = 26967) %>% 
  sf::st_buffer(dist = units::set_units(400, "m")) 

MARTA_buffer_group <- MARTA_buffer %>% 
  group_by(route_type) %>% 
  summarise()

pal <- colorFactor(palette = c("red", "yellow", "blue"), domain = MARTA_buffer_group$route_type)

leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addPolygons(data = MARTA_buffer_group %>% st_transform(crs = 4326), col = ~pal(route_type),
              popup = MARTA_buffer_group$route_type) 

Merge frequency by stop information with ACS data

stop_tract <- acs2020c %>% 
  st_transform(crs = 4326) %>% 
  st_join(stop_freq_sf) %>% 
  group_by(GEOID) %>% 
  summarise(n_departures = sum(n_departures),
            hhinc = mean(hhinc),
            r_tot = mean(r_tot),
            pct_wh = mean(pct_wh),
            pct_bl = mean(pct_bl),
            pct_novhc = mean(pct_novhc),
            r_tot = mean(r_tot),
            ln_pop_den = mean(ln_pop_den))

knitr::kable(head(stop_tract))

z <- function(vector){ round((vector - mean(vector))/sd(vector), 3) }

pass_per_stop_pal <- colorQuantile(palette = "Reds", domain = stop_tract$n_departures)

leaflet() %>% 
  addProviderTiles(providers$CartoDB.DarkMatter) %>% 
  addPolygons(data = stop_tract, 
              fillColor = pass_per_stop_pal(stop_tract$n_departures),
              fillOpacity = 0.8,
              color = 'grey',
              weight = 1, 
              popup = paste0("ID: ", stop_tract$GEOID, ", Value: ", stop_tract$n_departures)) %>% 
  addLegend(position = "bottomright",
            values = stop_tract$n_departures,
            pal = pass_per_stop_pal)

message(str_c("there are ", sum(is.na(stop_tract$n_departures)), " NAs in n_departures"))

## there are 84 NAs in n_departures

plotly::ggplotly(departure_count)